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Three terms, “Waterman’s T-matrix method", "extended boundary condition method 
(EBCM)”, and “null field method", have been interchangeable in the literature to indicate 
a method based on surface integral equations to calculate the T-matrix. Unlike the 
previous method, the invariant imbedding method (IIM) calculates the T-matrix by the 
use of a volume integral equation. In addition, the standard separation of variables 
method (SOV) can be applied to compute the T-matrix of a sphere centered at the origin 
of the coordinate system and having a maximal radius such that the sphere remains 
inscribed within a nonspherical particle. This study explores the feasibility of a numerical 
combination of the IIM and the SOV, hereafter referred to as the IIM + SOV method, for 
computing the single-scattering properties of nonspherical dielectric particles, which are, 
in general, inhomogeneous. The IIM + SOV method is shown to be capable of solving 
ligbt-scattering problems for large nonspherical particles where the standard EBCM fails 
to converge. The IIM + SOV method is flexible and applicable to inhomogeneous particles 
and aggregated nonspherical particles (overlapped circumscribed spheres) representing 
a challenge to the standard superposition T-matrix method. The IIM + SOV computational 
program, developed in this study, is validated against EBCM simulated spheroid and 
cylinder cases with excellent numerical agreement (up to four decimal places). In 
addition, solutions for cylinders with large aspect ratios, inhomogeneous particles, and 
two-particle systems are compared with results from discrete dipole approximation 
(DDA) computations, and comparisons with the improved geometric-optics method 
(IGOM) are found to be quite encouraging. 

© 2012 Elsevier Ltd. All rights reserved. 


1. Introduction 

In various scientific disciplines (bio-optics, photonics, 
astrophysics, and atmospheric radiative transfer and 
remote sensing), accurate and efficient computations of 
the optical properties of dielectric particles are often 
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required. The Lorenz-Mie theory and its modifications 
[1-5] applicable to homogenous or layered spheres cannot 
be used to compute the optical properties of morphologi- 
cally complex particulates. Great strides have been made 
toward accurate simulations of the scattering of light by 
particles of various shapes and/or chemical compositions, 
but, although a variety of methods have been developed 
[6,7], each method has its own strengths and weaknesses. 
For example, the full-electromagnetic wave methods of 
solving Maxwell’s equations become inefficient or even 
inapplicable when the particle size is excessively large, 
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whereas the semi-empirical geometric-optics method is 
applicable to large particles but fails for small particles 
when the “ray” concept is not valid. Tremendous effort has 
been expended to broaden the computational domain of 
the existing light scattering computational methods e.g„ 
[8-16] and to gain a better understanding of the single- 
scattering properties of nonspherical particles with size 
parameters ranging from the Rayleigh to geometric-optics 
regimes. 

Waterman’s T-matrix method (TMM) [17,18] is an 
accurate and powerful tool capable of yielding a highly 
accurate numerical solution for the scattering of light by 
nonspherical particles [19-22]. In the literature, the techni- 
que is sometimes referred to as either the extended bound- 
aiy condition method (EBCM) or the null field method. In 
contrast to many numerical techniques that explicitly con- 
sider the particle orientation and polarization state of the 
incident light in the simulations [23-28], the TMM calcu- 
lates the T-matrix, a quantity independent of the propaga- 
tion direction, incident light polarization state, and the 
scattering direction, and which allows for efficient compu- 
tation of the orientation-averaged optical properties [29]. 
The conceptual framework of the TMM has subsequently 
been expanded to handle composite particles, layered 
particles, and more complicated scattering cases e.g., 
[21,22,30-33]. These developments have demonstrated the 
EBCM to be just one possible path to compute the T-matrix 
of the scattering object. Some attempts [8,34-37] to 
improve the limited applicability of the EBCM to certain 
particles have focused on numerical instability, convergence 
issues, and loss of precision; however, the maximum size 
parameter value for a convergent EBCM solution strongly 
depends on the particle shape. Specific details of various 
TMM implementations and relevant applications can be 
found in the texts [20-22] and a reference database 
[38,39]. The cumulative body of relevant research has made 
the TMM one of the most widely used approaches to obtain 
highly accurate numerical optical properties of morpholo- 
gically complex particles with moderate aspect ratios and 
size parameters ranging from zero to ~ 200. 

As mentioned, the T-matrix, relating the incident to 
the scattered field expansions in vector spherical wave 
functions (VSWFs), can be computed from several alter- 
native approaches in addition to the EBCM. Johnson [40] 
derived the T-matrix from the standard electromagnetic 
volume integral equation (VIE) and developed an invar- 
iant imbedding method (IIM) to iteratively calculate the 
T-matrix. Schulz et al. [41] obtained the T-matrix of 
spheroids based on the separation of variables (SOV) 
method in spheroidal coordinates. A discrete dipole 
moment method to calculate the T-matrix was developed 
by Mackowski [42], while Loke et al. [43] incorporated the 
discrete dipole approximation (DDA) into the point- 
matching method to calculate the T-matrix. A superposi- 
tion T-matrix method (STMM) for multiple-sphere 
clusters developed by Peterson and Strom [30] and 
Mackowski and Mishchenko [44] is based on the addition 
translation theorem for vector spherical wave functions 
(VSWFs). In principle, any computational method that 
solves Maxwell’s equations can be employed to calculate 
the T-matrix, although the computational efficiency, the 


computer memory requirements, the complexity of the 
numerical implementation, and the range of practical 
applicability can be quite variable. 

The IIM for the calculation of the T-matrix has drawn 
scant attention since Johnson’s study [40] (with specific 
applications to relatively small particles) was published in 
1988. For example, the IIM is referenced neither in the 
T-matrix books [20-22] nor in the established T-matrix 
reference database [38,39]. Moreover, according to the ISI 
Web of Knowledge, Ref. [40] has previously been cited 
only 7 times. The use of the IIM in T-matrix calculations 
can be traced to its application to scattering problems in 
quantum mechanics [45]. Note that in addition to its 
application to the solution of Maxwell’s equations, the IIM 
has been applied to the solution of the radiative transfer 
equation [46[. During the 1990s and early 2000s, signifi- 
cant advancement in relevant numerical techniques and 
computer resources has been made, and we believe the 
time has come to revisit the IIM for state-of-the-art 
numerical implementations. 

This study explores the application of a numerical 
combination of Johnson’s IIM method and the SOV (here- 
after, IIM + SOV) to light scattering by large nonspherical 
and inhomogeneous particles. The remainder of the paper 
is organized as follows: Section 2 outlines the fundamen- 
tals of the T-matrix method and the invariant imbedding 
procedure; Section 3 includes the implementation of the 
IIM + SOV method to simulate light scattering by repre- 
sentative nonspherical and inhomogeneous particles, and 
both validates the accuracy and illustrates the efficiency 
by comparing IIM + SOV method results versus their 
counterparts computed from other methods; and. 
Section 4 summarizes our study. 


2. Theoretical basis 

2.1. The T-matrix 


To elaborate the concept of the IIM for the computation 
of the T-matrix, we begin with the definition of the T-matrix 
based on the expansion of the incident and scattered fields 
in terms of VSWFs. The definition of the T-matrix depends 
on the adopted functional basis. To facilitate a combination 
of the T-matrix computation and the analytical orientation- 
average algorithm outlined in Mishchenko et al. [21], we 
adopt the VSWFs in the exponential form in spherical 
coordinates. Let us consider the scattering of a plane 
electromagnetic wave E (r) by a finite volume with a 
non-unity relative refractive index surrounded by an infinite 
homogeneous, isotropic, and non-absorbing host medium. 
We expand the incident and scattered fields in terms of 
VSWFs as: 


OO n _ 

E'“(r)= ^ ^ YmnieMnir) 

n = 1 m = -n 


^mn 

^mn 


( 1 ) 


E^“(r)= £ ^ Y mnie, mn(r) 

n=im = -n 


Pmn 

Qmn 


r > r>. 


(2) 


where r> is the radius of the smallest circumscribed 
sphere of the scattering volume centered at the origin of 
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the coordinate system and the angular function Ymn(^.0) is 
a 3 X 3 matrix given by 




pinup 


■ 2n+\ 1’/^ 

4nn{n + 1) 

0 0 yn(n+T)do^(0)“ 

Tmn (0) 0 

~'i^mn{0) i'^mn(0) 0 


(3) 


where dg^{9) is the Wigner d function, while Ttmn and T^n 
are given, respectively, by 

Timn(d)= ^^[fig^om(0) = 2 \/n(n+ 1) [dj„(0) + d_ j^(6)] , (4) 


Timn(0) — ^dgujiG)— — \J n(n + 1) [dj^(d)— d_^^(6)] . (5) 

The radial functions J„(r) and H„(r) in Eqs. (1) and (2) 
are 3x2 matrices, given by 


J„ = 


Jn(fet) 

0 

0 


0 

YrWr [rjnikr)] 
y/n(n + 1)j„{kr)/kr 


(6a) 


H 


n — 


h0\kr) 

0 

0 


0 

[rh<’)(fcr)] 

^n(n+l)h6>(fer)/fcr 


(6b) 


where k is the wavenumber, jnikr) is the spherical Bessel 
function of the first kind, and hj’^kr) is the spherical Hankel 
function of the first kind. The following relations hold: 

Ymn Jn — YmnHn = {Mmn.Nmn)t (7) 


where Mmn and N^n are column vector spherical wave 
functions that satisfy the radiation condition at infinity, while 
RgMmn and RgNmn are those regular at the origin. Note that a 
bold symbol with one bar represents the column of {9,(p,r) 
components of a vector in the spherical coordinate system, 
while the symbol with double bars represent a matrix. 

The T-matrix is now defined by 


Pm'n' 

Qm'n' 


CO n 


E E 

n = 1 m = -n 


■ yl 1 
^ m'n'mn 


j21 

^ m'n'mn 


jl2 

^ m'n'mn 


j22 

^ m'n'mn 


^mn 

^mn 


( 8 ) 


Given the T-matrix, the amplitude scattering matrix 
and the phase matrix can be obtained. For a particle in a 
fixed-orientation with the incident light direction aligned 
with the z-axis, we need only to consider the T-matrices 
with m= + l. In this special orientation, for an axially 
symmetric particle, the amplitude scattering matrix is of 
the Lorenz-Mie type and given by 


where the two diagonal elements are 
°° 2n + 1 

Sz = E n{n+r) i'^m:n{cose) + bnTin{cose)], (10) 


Si = 


^ 2n + l 


[a„n„(cos6) + b„T„(cos6)] . 


( 11 ) 


In Eqs. (10) and (11), 6 is the scattering angle (defined 
as the angle between the incidence and scattering direc- 
tions), and the two coefficients are related to the T-matrix 
elements: 

n' = 1 


+ (13) 

If the incident light is along the negative z-axis direc- 
tion, we have 

- E + (14) 

n' = 1 


n' = 1 ’ 


For a particle that has mirror symmetry with respect to 
the xy-plane, Eqs. (14) and (15) are the same as Eqs. (12) 
and (13). In the simple case of a radially symmetric 
spherical particle, we have 


a„ = -T 


22 

Inin 


, fan = -T 


11 

Inin- 


(16) 


Note that the normalized phase matrix element P 22 /P 11 
derived from Eq. (9) is unity, even though the particle is 
nonspherical. For an axially symmetric particle with 
either a fixed or random orientation, the procedures 
necessary to calculate the amplitude scattering matrix 
and Mueller matrix from the T-matrix are described in 
detail by Mishchenko et al. [21]. In our study, we employ 
the analytical orientation average algorithm from 
Mishchenko et al. [21] and Mishchenko [29] in the 
computation of the phase matrix of randomly oriented 
particles, and the fixed-orientation algorithm from Mis- 
hchenko et al. [21] and Mishchenko [47] when the 
direction of incidence is not along the symmetry axis. 


2.2. The invariant imbedding method 

We recapture the basic principle of the IIM developed 
by Johnson [40] to calculate the T-matrix on the basis of 
an electromagnetic volume integral equation and sum- 
marize the major equations involved in the final compu- 
tation. An arbitrary scattering volume (homogenous or 
inhomogeneous) can be viewed as an inhomogeneous 
sphere (i.e., a portion of the sphere has the dielectric 
properties of the scattering particle and the reminder is 
medium/vacuum). In the spherical coordinate system, the 
particle can be discretized in terms of multiple inhomo- 
geneous spherical layers, as shown in Fig. 1(a). The 
principle of the IIM is to obtain the T-matrix of a larger 
sphere of p layers based on the T-matrix of a smaller 
sphere of p-1 layers. The initial T-matrix is zero at the 
origin of the coordinate system. 

To represent the T-matrix in a compact form, a combined 
index I is defined to represent two indices m and n via 
l=n(n+^)+m [48[. If the series expansions in Eqs. (1) and 
(2) are truncated to nmax. the range of 1 is [1, Imax] where 
lmax=nmaInmax+2). The T-matrix is now a l„,ax x Lax square 
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Fig. 1. Schematic geometry of an invariant imbedding approach to 
compute the T-matrix. (a) Discretize the nonspherical particles in terms 
of multiple-layered inhomogeneous spheres, (b) A combination of SOV 
and IIM to compute the T-matrix. (c) The physical interpolation of the 
IIM equation [17]. 

supermatrix whose element is a 2 x 2 matrix defined in Eq. 
(8). The fundamental equation to calculate the T-matrix, 
given by Johnson [40], is 

f(''p) = Qii(rp)+[I + Q,2(rp)][I-f(rp_])0;22(rp)]-^ 

><T(rp_i)[I + Q,2i (rp)], (17) 

where Pp is the radius of pth layer, 

=11 =12 
T (Pp) T (Pp) 

=21 =22 
T (Pp) T (Pp) 


is the T-matrix of a particle composed of p layers, I Is a 
2lmax X 2lmax Unit matrix, and the supermatrices Q,,j (l„,ax x 
Imax) are defined hy 


Qii ('■p) 

= ifeT(rp)5(tp)J(rp), 

(18) 

Ql2(tp) 

= ifeT(rp)5(tp)H(rp), 

(19) 

52,1 (tp) 

= ikH^ (rp)5(rp)J(Pp). 

(20) 

52,2 (tp) 

= ikH^ (rp)5(Pp)H(rp). 

(21) 


T(Pp) = 


where j (Pp) and H (Pp) are Imax x Imax diagonal sup«'matrices 
with each diagonal element defined by Eq. (6). Q.(Pp) and 
the relevant quantities are given by 


Q(Pp) =Wp 


[I-WpU(rp)g(Pp,rp)], 


(22) 




7^ 

( 6 , 0 ) 


X [m(p)2 - 1 ]Z(p)Y f ,„,„,)(6, 0), 


Z(p) = 


l/m^ 

0 

0 


0 0 
1 0 
0 1 


gn(rj') = < 


ikE„{r)j „(p'); r > r' 

|[H„(p)f n(r')+jn('-)H^(r')]; r = r' ■ 
ifcJ„(r)H’^„(r'); r < r' 


(23) 


(24) 


(25) 


In Eq. (22), U is a Imaxxlmax supermatrix with each 
element defined in Eq. (23), g is a tmaxxlmax diagonal 
supermatrix with each element defined in Eq. (25), and w„ 
is the weight at each discrete radius. In Eqs. (23) and (24), 
m is the refractive index. Eqs. (17)-(25) are similar to 
those in Section E in Johnson [40]; however, they differ 
because our formalism is based on exponential <p depen- 
dence and Johnson employed even and odd functions. 
An explicit derivation of Eq. (17) based on the exponential 
cp dependence Is given in Appendix. 

The direct application of Eq. (17) to numerical compu- 
tation Is Inefficient. To overcome the computational inef- 
ficiency, the SOV method can be employed to calculate the 
T-matrix of the sphere centered at the origin and inscribed 
within the scattering particle to the maximum extent. The 
llM is applied in the volume between the circumscribed 
sphere and the inscribed sphere centered in the coordinate 
systems as shown in Fig. 1(b). The Lorenz-Mie coefficients 
are related to the T-matrlx elements through Eq. (16); 
therefore, Eqs. (16)-(25) provide a closed set of mathe- 
matical equations to solve Maxwell’s equations. 

Johnson [40] did not discuss the physical interpolation 
of Eq. (17), which is shown in Fig. 1(c). The contribution of 
the surface (Pp) independent to the scattering is given by 
0 , 11 , the reflection of the incident source and the source 
from the ^rface Pp at the surface Pp_i is given by 
T(Pp-i)[I +0,2,1 (fp)]. and the contribution of the source 
fromjhe surface Pp__i ajter interaction with the surface Pp 
is Qi 2 (tp)T(tp-i)[I + 02,i (tp)]- Multiple interactions 
between two surfaces are taken into account through 
matrix inversion; 


[I-T(Pp_l)Q,2.2(''p)] ’ 


OO 


E[T(rp-,)02,2(rp)]'. 

k = 0 


(26) 


Note that the arrows in Fig. 1(c) explain the interaction 
relation involved In Eq. (17) and should not be construed 
as rays of light in geometric optics. 

The only shape-and-refractive-index-dependent part 
of the computational program is the U-matrix defined in 
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Eq. (23), which can be written by the following explicit 
expression, 


Based on Eqs. (30) and (31), 

rll t22 


,, = 0. when(-l)"+"' = -l, 


(32) 


U„ 


= /c^r2(-l) 


'.m + m' 


2n + l 2n' + l 


47tn(n+l) 


47tn'(n' + l) 


: /'Zn pTl 

/ d(j} d6sin6exp[-i{m-m')(j)][8{r,6,(l))-\)] 
Jo Jo 


^mn (0)7Tm'n'(^) + Tmn 
i \j^mn Tmn {0)Km'n'{d)\ 

0 


1 \j^mn Tmn (0)7ltn'n'{^)\ 

^mn Tmn 


^nn'(n + l)(n' + l)do^ (e)do^,(6)/in'^(r,e,4)) 

(27) 


In comparison with the EBCM, modifying the composi- 
tion of the scattering particle without modifying the 
system of equations is simpler and significantly reduces 
the complexity of numerical implementation. 

The symmetry of the T-matrix can be explicitly under- 
stood from the U-matrix. If the particle is of axially rotational 
symmetry, after integration in terms of (p, Eq. (27) simplifies 
to 


, = = 0 when i— 1)*^^” =1 

I' ' mnmn' ''viitii y if — i. 


(33) 


In the case of a homogenous sphere, employing the 
orthogonal relationships of 


r 


d6sin6[7Zjnn{G)Tmn'{^) TmnW)'^mn'{G)] — 0» 


(34) 


desinO[e(r,e)--l)] 
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^mn (6)71 
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mn' (0)] 

i \j^mn ( 6 )t 

mn' { 9 ) + Tmn (6)71 

mn' (6)] Ttmn (6)71 

mn' (6) + Tmn (6)t 

mn' (9) 


0 


0 


^nn'(n + l)(n' + l)dj„ {0)da„, (8) 
e.(r,8) 


(28) 


Note that [e(r,6)-l)] is a step function in terms of 6 
and is zero when 6 is outside the particle. If Eq. (28) is 
valid at each discrete radius, the T-matrix must have the 
following property 


Tmnm'n' — ^mm'T'mnmn' > 

where we denote 


(29) 
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d0sin6[Kmn{d)'llmn'{d) + Tmn{&)Tmn'iG)] — 2n -\- 1 
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d0sin8d'ij9)dZ(9) = 


dnn', 


gives the U-matrix in a diagonal form: 


*■ mnm'n' — 





■1 
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0 ■ 

■ ■T'l 1 

* mnm'n' 

^12 

^ mnm'n' 

Umnm'n' — ^ ^^(£~^)^mm''5nn' 

0 

1 

0 

't'21 

^ mnm'n' 

’T’22 

^ mnm'n' 


0 

0 

l/e_ 


If the axially symmetric particles have mirror symme- 
try with respect to the xy-plane, the U-matrix further 
simplifies to 


(35) 

(36) 

(37) 


Realizing the T-matrix symmetry in the case of an 
axially symmetric nonspherical particle, Eq. (17) can be 
separated into n^ax + 1 equations according to the index m: 


Umnm'n' --dmm'k [n^n' + lJ /o ^ d6sin6 [fi(r,6)- 1)] 

tnn' [^mn(6)7I mn' (9) + Tmn (6)t 

mn' (0)J ~^Cnn' [‘^mn(d)T mn' (6) + Tmn (9)71 

mn' (6)] 

iCnn' [7tmn(6)Tmn'(6)-l-Tmn(6)7rmn'(6)] Cnn' [7tmn(6)7tnin'(6)-l-Tmn(6)Tmn'(6)] 


c„„, ^nn'(n + l)(n' + l)dj„(9)dj;„(0) 
e(r,fl) 


(30) 


where 

Cm,' = [1 +(-!)"■""'], C„„, = [1 +(-1)"+"' + '] 


— _m - — m - — 

Tmnmn' (tp) = Q 1 1 (tp) + [I + Q. 12 (tp)][I— Tmnmn' (tp_] ] 

— m — - — m 
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(31) 


(38) 
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For each m, the size of the T-matrix is 
2(nmax-m+l) X 2(nmax-in + l). If the incident direction is 
along the axis, m=l is sufficient to obtain the solution. For 
an arbitrarily incident direction (particle orientation) or 
randomly oriented particle, m ranges from 0 to n^ax- In 
the case of a sphere, Eq. (17) becomes rimax separate 
equations in which the T-matrix is a 2 x 2 diagonal matrix. 

2.3. Numerical implementations 

The truncation number and the number of spherical 
layers to discretize the particle are two fundamental 
parameters to guarantee the convergence of the 
T-matrix. In the numerical calculations, the two computa- 
tional parameters are increased such that the solutions 
converge with acceptable precision. At different spherical 
surfaces, the T-matrix is truncated to a different number 
(i.e., the remaining elements in the T-matrix not involved 
in the computation are assumed to be zero) to overcome 
the overflow problem and increase the computational 
speed. The Gaussian quadrature is employed for the 
computation of the U-matrix (28), and the number of 
division points is increased so that the numerical error is 
smaller than a prescribed value (e.g., 10^’°). 

Clebsch-Gordan coefficients are indispensable in the 
analytical orientation average algorithm. The two- 
directional recursive algorithm [21] is employed for fast 
computation. The center value of the index ri of 
C ^nmni.i-m In Eq. (5.122), given in Appendix D in [21], 
is problematic when n' is larger than ~150. Instead, a 
computational procedure is implemented to search the 
center of the classical region. After the correction, the 
random-orientation averaging algorithm is tested for a 
sphere with a particle size parameter larger than 500 and 
gives the same results as its counterpart computed from 
the Lorenz-Mie theory. 

Computer memory requirements increase dramatically 
with respect to the particle size parameter. To reduce the 
memory burden that limits the largest applicable size 
parameter, we: (a) store non-zero T-matrix elements 
based on Eqs. (32) and (33) when the particle has a 
mirror symmetry with respect to the xy-plane; (b) reduce 
the three-dimensional D arrays in Eqs. (5.126)-(5.130) of 
Ref. [21 ] to two dimensional arrays by switching the loops 
in terms of the indices s and n; (c) store non-zero values of 
the B-array defined in Eq. (5.122) of Ref. [21]; and, 
(d) store the T-matrix elements and B-array in single- 
precision variables, f^^nn, proven to be zero when n-i is 
odd because of the symmetry of indices inherent in 
Clebsch-Gordan coefficients and the T-matrix. After the 
T-matrix elements for each index m are obtained in 
double precision, the single-precision and double- 
precision array to store T-matrix elements and B-array 
are found to give the same results within several decimal 
places. The final memory demand for a particle size 
parameter of 500, defined in terms of the radius of a 
circumscribed sphere, is less than 3 GB. 

To apply the IIM T-matrix method to axially symmetric 
particles with large size parameters, we have parallelized 
the Fortran 90 program based on the Message Passing 
Interface (MPI) technology. The calculations of sub 


T-matrices are distributed to different processors accord- 
ing to the index m. As the index m increases, the size of 
the T-matrix decreases. As a result, the computation time 
significantly decreases for large values of m. When the 
incident light direction is aligned with the symmetry axis, 
the algorithm does not need to be parallelized, because 
solving the equation for m = l is sufficient to obtain the 
scattering solution. 

Furthermore, reciprocity conditions and other rela- 
tionships [49] are employed to check the T-matrix and 
expansion coefficients of the phase matrix. All the equal- 
ities are in full agreement and the expansion coefficients 
of the phase matrix passed the check of the subroutine 
HOVENR in the Mishchenko’s EBCM program. For non- 
absorbing particles, the single-scattering albedo is found 
to be slightly greater than unity; however, the value is 
accurate to a significant number of places. 


3. Results and discussion 


3.1. Testing against EBCM: spheroids and cylinders 


We apply the IIM + SOV method to homogenous spher- 
oids and cylinders whose optical properties have been 
extensively studied by using the standard EBCM. A com- 
monly used EBCM code for spheroids and cylinders 
written by Mishchenko and Travis [50] is publicly avail- 
able at http://www.giss.nasa.gov/~crmim. The maximum 
size parameter handled by the EBCM decreases as the 
aspect ratio becomes extreme (i.e., very long or flat 
spheroids). In the extended-precision calculations, the 
maximum size parameter is larger than in the double- 
precision calculations. We have tested the IIM + SOV 
program (double-precision variables are used in the 
calculation of the T-matrix) in various cases calculated 
by the EBCM and excellent agreement between the two 
methods has been obtained. 




^ 0.5 
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Fig. 2. Comparison of phase matrix elements simulated from the EBCM 
and the IIM+SOV for randomly oriented spheroids. 
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Table 1 

The values of six nonzero phase matrix elements of randomly oriented prolate spheroids at 9 scattering angles simulated from the EBCM and the 
IIM+SOV. The "Diff’ rows in the table cells indicate the differences between the results computed from the two methods. 


en 

Pii 

P 22 

P 33 

P44 

P 12 

P 34 

0 

EBCM 

16.3398 

16.2871 

16.2871 

16.2345 

0.0000 

0.0000 

SOV + IIM 

16.3415 

16.2919 

16.2919 

16.2423 

0.0000 

0.0000 

Diff 

-0.0017 

-0.0048 

-0.0048 

-0.0078 

0.0000 

0.0000 

30 







EBCM 

4.9774 

4.9449 

4.9043 

4.8938 

-0.2142 

0.5491 

SOV+llM 

4.9772 

4.9452 

4.9045 

4.8943 

-0.2074 

0.5525 

Diff 

0.0002 

-0.0003 

-0.0002 

-0.0005 

-0.0068 

-0.0034 

60 







EBCM 

0.2987 

0.2835 

0.2156 

0.2257 

0.0976 

0.0596 

SOV+llM 

0.2989 

0.2836 

0.2151 

0.2252 

0.0958 

0.0590 

Diff 

-0.0002 

-0.0001 

0.0005 

0.0005 

0.0018 

0.0006 

90 







EBCM 

0.1459 

0.1220 

0.0755 

0.0970 

-0.0471 

-0.0340 

SOV + IIM 

0.1457 

0.1217 

0.0759 

0.0972 

-0.0470 

-0.0340 

Diff 

0.0002 

0.0003 

-0.0004 

-0.0002 

-0.0001 

0.0000 

120 







EBCM 

0.0621 

0.0356 

0.0027 

0.0272 

-0.0094 

0.0103 

SOV + IIM 

0.0620 

0.0355 

0.0025 

0.0269 

-0.0094 

0.0103 

Diff 

0.0001 

0.0001 

0.0002 

0.0003 

0.0000 

0.0000 

150 







EBCM 

0.0326 

0.0244 

-0.0197 

-0.0132 

0.0024 

-0.0058 

SOV + IIM 

0.0328 

0.0244 

-0.0198 

-0.0132 

0.0024 

-0.0058 

Diff 

-0.0002 

0.0000 

0.0001 

0.0000 

0.0000 

0.0000 


180 

EBCM 

0.0585 

0.0329 

-0.0329 

-0.0074 

0.0000 

0.0000 

SOV+IIM 

0.0583 

0.0327 

-0.0327 

-0.0071 

0.0000 

0.0000 

Diff 

0.0002 

-0.0002 

-0.0002 

-0.0003 

0.0000 

0.0000 


Table 2 

The extinction efficiency, the scattering efficiency, the absorption effi- 
ciency, the single-scattering albedo, and the asymmetry factor simulated 
from the EBCM and the IIM-tSOV. The "Diff” row in the table cells 
indicates the differences between the results computed from the two 
methods. 



Qext 

Qsca 

Qabs 

CO 

< cos 6 > 

EBCM 

3.2854 

2.2903 

0.9951 

0.6971 

0.8165 

SOV+IIM 

3.2850 

2.2901 

0.9949 

0.6971 

0.8165 

Diff 

0.0004 

0.0002 

0.0002 

0.0000 

0.0000 


Table 1 presents the comparison of the six nonzero phase 
matrix elements of a prolate spheroid in random orientation 
computed from the IIM+SOV (the second row for each 
scattering angle) and the benchmark results (the first row 
for each scattering angle) computed from the EBCM [29]. The 
ratio of semi-major axis (b) over semi-minor axis (a) is 2, the 
size parameter (kb) is 5.5, and the refractive index is 
1.5 + 0.11. The IIM+SOV can reproduce up to four decimal 
places depending on the scattering angle. In the computation, 
a Romberg integration technique is employed to improve the 
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Fig. 3. Comparison of phase matrix elements simulated from the EBCM 
and the IIM+SOV for randomly oriented compact cylinders. 
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Table 3 

Similar to Table 1, but for randomly oriented cylinders. 


en 

Pu 

P 22 

P 33 

P44 

P 12 

P 34 

0 

EBCM 

9.8167 

9.7796 

9.7796 

9.7426 

0.0000 

0.0000 

SOV-hI!M 

9.8176 

9.7817 

9.7817 

9.7458 

0.0000 

0.0000 

Diff 

-0.0009 

-0.0021 

-0.0021 

-0.0032 

0.0000 

0.0000 

30 







EBCM 

4.5720 

4.5226 

4.5064 

4.5045 

-0.2996 

0.2008 

SOV -h IIM 

4.5732 

4.5242 

4.5085 

4.5068 

-0.2932 

0.2010 

Diff 

-0.0012 

-0.0016 

-0.0021 

-0.0023 

-0.0064 

-0.0002 

60 







EBCM 

0.8046 

0.7415 

0.7084 

0.7536 

-0.0769 

0.1585 

SOV -h IIM 

0.8043 

0.7411 

0.7075 

0.7532 

-0.0765 

0.1581 

Diff 

0.0003 

0.0004 

0.0009 

0.0004 

-0.0004 

0.0004 

90 







EBCM 

0.1972 

0.1442 

0.1151 

0.1595 

0.0367 

-0.0007 

SOV -h IIM 

0.1973 

0.1440 

0.1148 

0.1594 

0.0363 

-0.0009 

Diff 

-0.0001 

0.0002 

0.0003 

0.0001 

0.0004 

-0.0002 

120 







EBCM 

0.1080 

0.0840 

-0.0193 

-0.0006 

-0.0106 

-0.0479 

SOV-hllM 

0.1077 

0.0835 

-0.0190 

-0.0002 

-0.0108 

-0.0481 

Diff 

0.0003 

0.0005 

-0.0003 

-0.0004 

0.0002 

0.0002 

150 







EBCM 

0.1226 

0.1061 

-0.0964 

-0.0819 

-0.0042 

-0.0146 

SOV-hllM 

0.1223 

0.1057 

-0.0960 

-0.0815 

-0.0043 

-0.0145 

Diff 

0.0003 

0.0004 

-0.0004 

-0.0004 

0.0001 

-0.0001 

180 







EBCM 

0.1731 

0.1261 

-0.1261 

-0.0790 

0.0000 

0.0000 

SOV -h IIM 

0.1730 

0.1259 

-0.1259 

-0.0788 

0.0000 

0.0000 

Diff 

0.0001 

0.0002 

-0.0002 

-0.0002 

0.0000 

0.0000 


Table 4 

Similar to Table 2, but for randomly oriented cylinders. 



Qext 

Qsca 

Qabs 

CO 

<cos 0} 

EBCM 

2.5322 

2.4471 

0.0851 

0.9664 

0.7090 

SOV + IIM 

2.5313 

2.4462 

0.0851 

0.9664 

0.7092 

Diff 

0.0009 

0.0009 

0.0000 

0.0000 

-0.0002 


numerical accuracy of the T-matrix until the decimal points 
in the table will not change by increasing the number of 
spherical layers to discretize the spheroid. The truncation 
number is increased in order for the final results to converge. 
Although some decimal points vary between the results 
computed from the two different methods, the forward and 
backscattering phase matrix computed from the two meth- 
ods satisfies the inherent relations: — ^11+^22 + ^33 — ^44=0 
at the fon/vard scattering angle and Pii-P 22 +f’ 33 -f’ 44=0 at 
the backscattering angle [21]. Table 2 shows a comparison 
between the extinction efficiency, the scattering efficiency, 
the absorption efficiency, the single-scattering albedo, and 
the asymmetry factor computed by the two methods. The 
results agree up to four decimal places. 

Fig. 2 shows a comparison between the phase matrix 
of a large randomly oriented prolate spheroid computed 


with the EBCM and with a combination of IIM and SOV. 
The ratio of the semi-major to the semi-minor axis is 2. 
The size parameter, defined in terms of the semi-major 
axis, is 95. The refractive index is chosen to be 1.311, 
which was used to test the computational capabilities of 
the EBCM (e.g., [34,50]) and is used as a reference in the 
present study. The comparison exhibits virtually the same 
results. The computational time is approximately 22 min 
using 32 (2.8 GHz) processors for the IIM + SOV method 
and 13 min using a single (2.8 GHz) processor for the 
extended-precision EBCM. Thus, using the EBCM is highly 
recommended within its applicable size parameter 
regime due to its significantly shorter processing time. 

Numerical simulations with cylinders have also been 
carried out, similar to those for spheroids. Tables 3 and 4 
show the optical properties of a randomly oriented circular 
cylinder with an equal-surface-area-sphere size parameter of 
3, a diameter-to-height ration of 0.5, and a refractive index of 
1.53 + 0.008i [51[. Similar accuracy to that illustrated in 
Tables 1 and 2 is obtained. Fig. 3 shows the comparison of 
six independent non-zero phase matrix elements for a large 
randomly oriented cylinder simulated from the EBCM and 
the IIM+SOV for a size parameter of {kD/2 = 50, kH 12 = 50}, 
where D is the diameter and H is the height. As demonstrated 
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in Fig. 3, an overall agreement is achieved but small 
differences are noticeable in the P43 element. A halo near 
46°, commensurate with geometric optics, is evident at this 
size parameter. A discussion of the particle size necessary to 
produce halos has been reported where the phase function of 
a cylindrical particle was computed with the EBCM for size 
parameters of 40, 80,120, and 180, with the size parameter 
defined in terms of kD/2 [52], 


3.2. Comparison with the improved geometric optics 
method: large spheroids and cylinders 


Fig. 4 shows the application of the IIM + SOV method to 
a size parameter of 300, i.e., where no other rigorous or 

10' 

10’ 

10“ 

10= 

Q.’” 

10' 

10 ° 

10 "’ 
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numerical methods are currently applicable. The simu- 
lated results are compared with the approximate 
geometric-optics results computed with the Improved 
Geometric-Optics Method (IGOM) [14]. As demonstrated 
in Fig. 4, the geometric-optics results are similar to those 
computed from the IIM + SOV. Note that the interference 
among rays is neglected in the IGOM calculations and 
explains the missing oscillations in the IGOM results. 
According to Mishchenko and Travis [50], for the 
extended-precision EBCM, the maximum convergent size 
parameter of a prolate spheroid with an aspect ratio of 2 
is 112 for the semi-major axis. Similar to Fig. 4, Fig. 5 
shows an oblate spheroid case. In Ref. [34], the EBCM was 
applied to large oblate spheroids with semi-major axis 
size parameters up to 102.3216, which has a surface- 
equivalent-sphere size parameter of 85. The present size 
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Fig. 4. Comparison of phase matrix elements simulated from the Fig. 6. Comparison of phase matrix elements simulated from the 

IIM+SOV and the IGOM for randomly oriented spheroids. IIM+SOV and the IGOM for randomly oriented prolate cylinders. 
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Fig. 5. Similar to Fig. 4, but for oblate spheroids. 


Fig. 7. Similar to Fig. 6, but for oblate cylinders. 
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parameter of the spheroid is larger by a factor of three 
than that handled by the EBCM. To obtain a better 
knowledge of the accuracy of the geometric-optics 
method for smaller size parameters, the reader is referred 
to the comparison of the complete phase matrix elements 
of a randomly oriented spheroid simulated from the 
conventional geometric optics method (CGOM) and the 
EBCM in [34,53], or the IGOM and the EBCM in [14], 

The comparisons of the phase matrix elements simu- 
lated from the IIM+SOV and the IGOM of prolate and 
oblate cylinders are presented in Figs. 6 and 7, respec- 
tively. The long-to-short-dimension ratio is 2, and the size 
parameter defined in terms of half of the long dimension 
is 300. The maximum convergent size parameter for an 
oblate cylinder of the extended-precision EBCM at this 
aspect ratio is 70 [50]. The excellent agreement between 
the IIM+SOV and the IGOM suggests the high reliability 
of the geometric optics method for large size parameter 
ranges. The differences between the IGOM and IIM+SOV 
for spheroids are relatively larger than those for cylinders. 
One possible explanation is the local radius of curvature 
on a cylindrical surface is much larger than on spheroidal 


surface and better fits the conditions of validity of the 
eikonal approximation for geometric optics. For cylinders, 
the radii of curvature on the top and bottom faces and the 
curves on side faces along the height are infinite, and the 
radius of curvature of the circular cross section is propor- 
tional to the size parameter. However, for spheroids, the 
radius of curvature is largest at the center and decreases 
as either end is approached. 

3.3. Comparison with discrete-dipole-approximation 
method: extreme cylinders, inhomogeneous particles, and 
aggregates 

We have shown the IIM+SOV method to be applicable 
to high aspect ratio and large size parameter cylindrical 
particles, which are out of the current computational range 
of the EBCM. For the purpose of comparison, we simulated 
the same scattering cases by using the DDA method on 
parallel computer clusters. To facilitate the DDA simula- 
tion for large size parameters, we considered a fixed 
orientation and the direction of the incident light aligned 
with the symmetry axis. Fig. 8 shows the comparison of 
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Fig. 8. P,i elements computed from the ADDA and the IIM+SOV for cylinders of large aspect ratios where the standard EBCM fails. 
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Fig. 9. Phase functions for six types of inhomogeneous particles. 


the phase-matrix elements for circular cylinders of various 
aspect ratios simulated from the ADDA and the IlM+SOV. 
The fundamentally different algorithms give almost the 
same results. The advantage of the T-matrlx over the ADDA 
is the computational efficiency of the llM + SOV in cases of 
random orientations and large size parameters. 

The computation of the T-matrix from the llM+SOV 
has been demonstrated to be robust in handling light 
scattering by nonspherical particles. In addition to non- 
spbericity, natural particles often have heterogeneous 
compositions. For example, Aslan dust aerosols originating 
in desert source regions can be coated by sulfates or soot 
when passing through heavily polluted downwind indus- 
trial regions [54], and Saharan dust grains transported over 
the Atlantic Ocean are often coated by sea salt. 

As Illustrative examples, we considered sk types of 
representative inhomogeneous particles whose phase func- 
tions are shown in Fig. 9. In the first column, the size 
parameters of a sphere (defined in terms of its radius) and 
a prolate spheroid of aspect ratio 2 (defined in terms of its 
semi-major axis) are 30. The refractive indices of the lower 
and upper part are 1.33 and 1.53. In the second column, the 
size parameter of a sphere (defined in terms of radius) is 15 
and a prolate spheroid of aspect ratio 2 (defined in terms of 
its semi-major axis) is 30. The ratio of the size parameters of 
the outer and inner radii is 0.8. The azimuthal angles 
necessary to separate different compositions are 30° and 
150°. The core refractive index is 1.53, the black part Is 1.33, 
and the gray part Is 1.70. Except for the refractive index, the 
size parameter and geometries in the third column are the 


same as their counterparts in the second column. The size 
parameter of the spherical inclusion is 2.5 in the sphere, and 
7.5 in the spheroid. The distance between the center of 
inclusion and the host particle is 10. In the simulations, the 
refractive index of the host particle is 1.33 and of the 
inclusion is 1.53. In the ADDA simulations, the more inho- 
mogeneous the particle the slower the convergence, and to 
guarantee numerical accuracy, the number of dipoles per 
wavelength needs to be larger. Note that light scattering by a 
sphere with an inclusion was solved by Vldeen et al. [55] 
using an extension of Lorenz-Mie theory. The advantage of 
the IIM+SOV method is the simplicity with which it deals 
with diverse inhomogeneous particles. 

Natural particulate matter has a tendency to aggregate 
and form single nonspherical particles. Scattering calcula- 
tions for multiple nonspherical particles have been 
reported [56-59] based on the STMM. Fig. 10 shows 
phase functions of several two-particle systems computed 
by the IIM + SOV method and the ADDA. The shaded 
component represents a different refractive index from 
the other component. In the first column, we use the IIM 
to compute the T-matrix, although the IIM + SOV can be 
employed if the center of the particle geometry is shifted. 
In the remaining cases, we use the IIM+SOV because the 
inscribed sphere is easily Identified. Note that for the bi- 
sphere and bi-prolate-spheroid systems, it is possible to 
compute the T-matrix for individual components and to 
use the STMM to obtain the T-matrix of the cluster. 
However, for the case shown in the third column, the 
superposition STMM falls because the circumscribed 
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Fig. 10. Phase functions computed from the IIM+SOV and the ADDA for geometries indicated within the figure. 


spheres of the two spheroids overlap [21 ]. To compute the 
phase function as well as for simplicity, we consider the 
direction of the incident light to be aligned with the 
symmetry axis. A very good agreement is demonstrated. 
The size parameters are indicated in the upper panel. The 
aspect ratios for the spheroids (major axis/minor axis) are 
assumed to be 2, and the refractive index is 1.33. The 
geometry and size of the particles in the lower panel are 
the same as those in the upper panels except that the 
refractive index in the shaded areas is 1.53. 

4. Concluding remarlts 

The IIM + SOV method provides an alternative approach 
to solving light scattering by nonspherical inhomogeneous 
particles. As compared with the EBCM, Johnson’s llM 
method had received less attention, and its modeling 
capabilities had not been extensively explored. In addition 
to the theoretical formulation, we have provided physical 
interpretation of terms in the iterative T-matrix equation. In 
this study, the IIM+SOV method has been applied to 
challenging or unprecedented computational domains, such 
as large size parameters, extreme geometries, and densely 
packed nonspherical particles. 

The comparisons of the results simulated from the 
IIM + SOV, the EBCM, and the ADDA methods have illu- 
strated that the present numerical implementation of the 
IIM + SOV is correct. The simulated results satisfy the 


inherent relations, summarized by van Der Mee and 
Hovenier [49], with acceptable accuracy. However, a more 
accurate quantification of the applicability domain of the 
method for miscellaneous particles is still necessary. For 
example, for size parameters into the thousands, the 
algorithm must be parallelized according to a shared 
memory model. The extension of the IIM + SOV method 
to particles that lack axial symmetry is may also be a 
future research topic. 

Compared to the EBCM, the particle morphology in the 
computational program is profoundly easier to modify, 
especially in the case of layered and composite particles. 
To further enhance the efficiency and the numerical 
capabilities of the IIM + SOV T-matrix method, an integra- 
tion of the IIM, the EBCM, the SOV, and the superposition 
T-matrix method has been developed, and relevant results 
will be published in a separate paper. 
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Appendix A. Iterative T-matrix equations for 
inhomogeneous spherical layers 


Substituting Eqs. (46) and (25) into Eq. (44), we obtain 

n = im = -n 

(47) 


The electromagnetic volume integral equation as given 
in Refs. [7,21,60], reads 


£(r) = (f) + k^ j (m^-1)G(r-r')E(r)d^r. (39) 

where £'"‘^is the incident electric field and G(r-r ) is the 
dyadic Green’s function [48]: 

exp I ik|r-r' | 

G=tf+\v v). — V ( 40 ) 

^ ^ 47t|r-d| 

Let |l^-r'|=R; then the explicit expression for G is 
given by 


exp(ikR) 
^ 47tR 



1-ikR Ij 3RR 
(kRf \ ~ (kRf 


(41) 


and in a simplified form 

tm'n'(r) = Ym.„.(0,0)J„,(r)+ f dr' ^ ^ ym’nf0.ll>)g„(r.R)f^„m„fr’) 
■'0 = 

(48) 

fmnm'n-ir) = [ dQ'7 {O' ,<()') u(r)Z(r (f) (49) 

Jv 

Substituting Eq. (48) into Eq. (49) yields 
Fmnm'n'(iO “ ^mnm'n' ir)in’(r) 

pR oo n 

+ „ E U/nn.m n(r)gfi(r,r')F 

mnm'n' if) 


n = 1 rh = -n 


(50) 


where 


Umnm'n'(t) — / dl2[Ymn(d,(/>)]^’*'u(r)Z(r)Ym'n'((('0)- (fif) 


The singularity term associated with the 1/R^ depen- 
dence in Eq. (41 ) is as follows. 


*-*S 

G = 


1 

471 



1 

47tR 


-y-kR-^ 


]GR^ 


rR^j 


(42) 


where the cap denotes the unit vector in a specific 
direction. To treat the singularity, Eq. (39) can be written 
as 


E{r) = Ei„c(r) + k^ 



G(f-r') 


E(r')d^r\ (43) 

where G(r-r ) is the Green function valid at the source 
free regions. Eq. (43) is equivalent to 


+ |y<5(r-r ) 
K 


E'^ff(r) = E'"‘^(f)+ f dVG(r',r')[m(r')2-l]Z(r')£"-'^(r') (44) 

Jv 

whpre (r ) is the total field in source free regions, and 
Z(r') [40] is a Cartesian tensor given by 


Z(r )=J/tfi^fr+99+(f>^. (45) 

In the present context of a spherical coordinate sys- 
tem, it is preferable to classify the light scattering by 
particles into two groups: homogenous spheres and 
inhomogeneous spheres. The scattering by nonspherical 
particles is equivalent to that by an inhomogeneous 
circumscribed sphere (the refractive index of the space 
between the spherical surface and the particle space is 
different from that inside the particle). On the basis of 
Eq. (7), the dyadic Green’s function in the matrix form is 
written as 

^ CO n ^ j 

Co(f,r')= ^ ^ 7 mn{9,(l>)J^(f,r )\ ^„{9' ,4>'), (46) 

n = 1 m = -n 

where g„(f,r ) is a 3-by-3 matrix given by Eq. (25). 


Based on Eq. (48), the T-matrix is written as, 

= ^ fR = ^ 

E mnm'n' (T) = ik / ddj^ (t')E n2nm'n'(T')f 


(52) 


To numerically compute the T-matrix, the integral is 
represented as summations 


T(r„) = ik^w/(rj)F(n|rj), 


7 = 1 


(53) 


F(n|ri) = U(r,)J(ri)-l- ^ WjU(r,)g(ri,rj)F(n|rj). 

7 = 1 

Based on Eq. (54), we have 

_ _ n-1 _ 

F(n|r„) = U(r„)j(r„)+ ^ WjU(r„)g(r„,rj)F (n|rj) 
7 = 1 

+ Wn U(r„)g(r„,r„)F(n|r„), 

and 


l-w„U(r„)g(r„,r„)]F ( n|r„) = U(r„) 

n-1 


(54) 


J(tn)+ E'’^7^(''"’'7')f("l'7') 

7 = 1 


Inserting Eq. (25) into Eq. (56) gives 

F(n|r„) = w„-'w„ 

[/-w„U(r„)g(r„,rn)] 

_ n-1 

J(r„) + H(r„)ik ^ WjJ(rj)F{n\rj) 
7 = 1 

Define 


(55) 


(56) 


(57) 


Q = -^ q=ikJ2wjf {rj)F{n\rj), (58) 

[I-w„U(r„)g(r„,r„)] 
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then 

F(n|r„) = w„-’Q;|J(r„)+H(r„)g]. 

Based on Eqs. (53), (58) and (59), we obtain 

— n-1 = _ — 

T(r„) = ik Y) Wjf (r,)F(n|r,) +!kwJ^(rn)F(n|r„) 

j = i 

= f+ifcwJ^(r„)F(n|r„) 


(59) 


= q + ikf(r„)(l 


j(r„) + H(r„)q 


— 9+0-11 +Oi29 


— Q.11 +(1+0.12)9’ 

(60) 

[10] 

where Q, , = I'feTO J and Q,, 2 = ikf Q H. Beginning with 

[11] 

the following equation. 


T(r„.^i) = ik Y Wjf {rj)f{n--l \rj), 
j = i 

(61) 

[12| 

let 



f(n\rj) =f{n--l\rj)d+n 

(62) 

[13| 

then 


[14] 

T(r„_i)(I+p) = 9. 

(63) 


Based on Eqs. (54) and (62), we have 


[15| 

_ _ _ n-1 

F(n-1 |r,)(I +p) = D(ri)J(n)+ Y w,U(r,)f (n.rj) 
j=l 


[16| 

xF(n-l rj)(I+p)+WnD(r,)f(ri,r„)F(n|rn) 

(64) 

[17] 

Note that Eq. (54) is satisfied if n is replaced with n-1. 

[18] 

and Eq. (64) is in the following form 


[19] 

U(n)J(n)F = U(r,)f(ri,r„)0 [j(r„) + H(r„)q]. 

(65) 

[20] 

Comparing Eqs. (25) and (65), we obtain 


[21] 

P = 02,1 +02,29 

(66) 

[22] 

Combining Eqs. (63) and (66), we have 



_ _ 


[23] 

T(rn-l)(l + 02,l +02,29) = 9 

(67) 


where Q .2 1 = ifeH^OJ. and 0 , 2,2 = OH. 
Eq. (67) in terms of q, 

Solving 

[24] 

9 = [I~T(r„_i)02,2] 'T(rn_i)(I + Q2,i). 

(68) 

[25] 

Substituting Eq. (68) into Eq. (60), we have Eq. (17), 

[26] 

which is the final iterative formula to calculate the 

T-matrix . The information about particle shape and 

[27] 

composition is contained in the U-matrix given by Eq. 
(23). Note that the definition of the U-matrix (relevant 
quantities) is different if VSWFs are formulated based on 

[28] 

odd and even modes [40], although the same 

symbols 

[29] 

are used. 


[30] 
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